Se pretende estimar \(\mathbb{E}_{f}[g(X_{1}, X_{2})]\), donde \[ f(x_{1}, x_{2}) \propto \begin{cases} \mathrm{e}^{-4 (x_{2} - x_{1}^{2})^{2} + (x_{2} - 1)^{2}} &\text{si } x_{2} \leq 2, \\ 0 &\text{en otro caso.} \end{cases} \qquad g(x_{1}, x_2) = \sqrt{x_{1}^{2} + x_{2}^{2}} \]

library(plotly)

f_u <- function(x1, x2) {
  exp(-4 * (x2 - x1^2)^2 + (x2 - 1)^2) * (x2 <= 2)
}

malla <- data.frame(
  x1 = (-250:250) / 100,
  x2 = (-250:250) / 100
)
malla |>
  plot_ly(x = ~x2, y = ~x1) |>
  add_surface(z = outer(malla$x1, malla$x2, f_u))

Vamos a realizar la estimación requerida a partir de un paseo aleatorio de Metropolis, haciendo uso del paquete mcmc. Para ello, en primer lugar debemos definir el logaritmo de la densidad objetivo como una función que toma como argumento un vector con los dos valores \(x_{1}\) y \(x_{2}\).

log_f_u <- function(x) {
  x1 <- x[[1]]
  x2 <- x[[2]]
  if (x2 <= 2) {
    -4 * (x2 - x1^2)^2 + (x2 - 1)^2
  } else {
    -Inf
  }
}

A continuación, iniciamos el paseo aleatorio de Metropolis a partir del estado \(x_{1} = x_{2} = 0\), que según el gráfico se encuentra cerca de donde la densidad toma su máximo. Podemos comprobar que la función metrop devuelve la realización de la cadena de Markov en forma de matriz con dos columnas, ya que cada estado es una tupla de dos números reales.

library(mcmc)

paseo_Metropolis <- metrop(log_f_u, initial = c(0, 0), nbatch = 1e4)
head(paseo_Metropolis$batch)
##           [,1]       [,2]
## [1,] 0.0000000  0.0000000
## [2,] 0.0000000  0.0000000
## [3,] 0.0000000  0.0000000
## [4,] 0.7588121 -0.2884623
## [5,] 0.5888510 -0.2832625
## [6,] 0.5888510 -0.2832625

La regla empírica establece que para acelerar la convergencia de la cadena de Markov a la densidad objetivo, el porcentaje de estados aceptados debería ser del 25 %.

paseo_Metropolis$accept
## [1] 0.2293
paseo_Metropolis <- metrop(paseo_Metropolis, scale = .9)
paseo_Metropolis$accept
## [1] 0.2662
paseo_Metropolis <- metrop(paseo_Metropolis, scale = .95)
paseo_Metropolis$accept
## [1] 0.2416
paseo_Metropolis <- metrop(paseo_Metropolis, scale = .93)
paseo_Metropolis$accept
## [1] 0.2471

Como diagnóstico de convergencia, los siguientes gráficos parecen indicar que la cadena de Markov ha convergido a la densidad deseada.

library(coda)

paseo_Metropolis$batch |>
  mcmc() |>
  traceplot()

library(ggplot2)

ggplot(
  as.data.frame(paseo_Metropolis$batch),
  aes(x = V1, y = V2)
) +
  geom_contour(
    data = expand.grid(x1 = (-250:250) / 100, x2 = (-250:250) / 100),
    mapping = aes(x = x1, y = x2, z = f_u(x1, x2))
  ) +
  geom_point(size = .5)

El gráfico de autocorrelación parece indicar que, al aplicar el método de las medias por lotes, la longitud de estos debería ser de al menos 50 estados.

g <- function(x) {
  x1 <- x[[1]]
  x2 <- x[[2]]
  sqrt(x1^2 + x2^2)
}

library(purrr)

paseo_Metropolis$batch |>
  array_branch(margin = 1) |>
  map_dbl(g) |>
  mcmc() |>
  autocorr.plot(lag.max = 100)

Estamos ya en condiciones de obtener la estimación requerida, así como un intervalo de confianza para la misma.

paseo_Metropolis <-
  metrop(paseo_Metropolis, nbatch = 1e3, blen = 1e2, outfun = g)

paseo_Metropolis$batch |>
  mcmc() |>
  autocorr.plot()

estimacion <- paseo_Metropolis$batch |>
  array_branch(margin = 2) |>
  map_dbl(mean)
error_estandar <- paseo_Metropolis$batch |>
  array_branch(margin = 2) |>
  map_dbl(var) |>
  magrittr::divide_by(paseo_Metropolis$nbatch) |>
  sqrt()

probabilidad_cobertura <- 0.95
alfa <- 1 - probabilidad_cobertura
percentil <- qnorm(1 - alfa / 2)
intervalo_confianza <- estimacion + c(-1, 1) * error_estandar * percentil

Se obtiene entonces una estimación de 0.8557277, con (0.8365232, 0.8749322) un intervalo de confianza con probabilidad de cobertura 0.95.